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We study numerically the zero temperature Random Field Ising Model on cubic lattices of various 
linear sizes L in three dimensions. For each random field configuration we vary the ferromagnetic 
coupling strength J. We find that in the infinite volume limit the magnetization is discontinuous in 
J. The energy and its first J derivative are continuous. The approch to the thermodynamic limit 
is slow, behaving like L~ p with p ~ .8 for the gaussian distribution of the random field. We also 
study the bimodal distribution hi = zth, and we find similar results for the magnetization but with 
a different value of the exponent p ~ .6. This raises the question of the validity of universality for 
the random field problem. 
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^j- , I. INTRODUCTION 

o 

Despite its long history, the Random Field Ising Model (RFIM) is not yet fully undesrstood Q. It is known that 
— ^, \ in three dimensions and for weak disorder there is a paramagnetic to ferromagnetic phase transition but the nature 
of the transition is still unknown. What is most remarquable is that the RFIM is (together with branched polymers 
H) one of the very few cases where Perturbative Renormalization Group (PRG) can be analyzed to all orders of 
1 , pertubation theory J3|,[I| and this analysis leads to wrong conclusions (while for the branched polymers the same 
analysis is correct!). It has been argued that this failure of PRG is due to replica symmetry breaking This 
failure raises the more general question of the validity of perturbative renormalization group for disordered systems. 

In the present paper we study numerically the phase transition of the RFIM in three dimensions. Several numerical 
studies have already been performed p|-p^| but the problem has been proven to be very difficult and, in our opinion, 
^ • no definite conclusion about the nature of the transition has been drawn yet. We simulated much larger sizes than 
before (up to 90 3 ) with much higher statistics and we establish new results on the phase transition of the RFIM. 
We find that the magnetization is discontinuous at the transition both for the gaussian and bimodal distribution. 
The energy and its derivative are continuous only for the gaussian distribution. The exponents are found to be non 
classical. Such a behaviour has already been suggested on the grounds of real space renormalization group 
(This result has been often overlooked in the past due to the approximations of the method.) We also find different 
exponents for different random field distributions, which is in contradiction with the PRG. 
The Hamiltonian of the RFIM is given by 

H = —J <7i<Tj — hi<Ji (1) 

<i,j> i 

where Yl^j runs over neighbouring sites of the lattice (we have only considered three dimensional cubic lattices with 
periodic boundary conditions) and {hi} are independent random variables identically distributed with zero mean and 
variance one. For a given random field sample, one can vary both the ferromagnetic coupling J and the temperature 
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T, i.e. the phase boundary is a line in the J, T plane. It is thought, in accordance with PRG, that the nature of the 
transition and the value of the exponents do not depend on the position on the transition line, nor on the direction 
one crosses it, and that this is true down to zero temperature. So it is advantageous to work at T = where it has 
been shown that the RFIM is equivalent to the problem of finding a maximum flow in a graph [fl5f , for which very fast 
(polynomial) algorithms are known. Note that these algorithms provide the exact groud-states and therefore there 
is no thermalization problem. Simulations using such algorithms have already been performed in the past |^o| , |ll| . 
In the present paper we use the latest version of the algorithm developped by Goldberg and Tarjan Jlq] , which we 
optimized for the case of the cubic lattice. It has been shown [[l6) that this algorithm converges to the ground state 
in a time t < L 6 lni where L is the linear size of the cubic lattice. We found experimentally that t ~ L A ]l7j ]. 



II. GAUSSIAN DISTRIBUTION OF RANDOM FIELD. 



We first consider the case of a gaussian distribution of the random field with variance equal to one. It is customary 
to present the data as a function of x = 1/J. For every sample {hi}, we chose equally spaced x's with 6x — .0125 and 
found the corresponding ground state. We considered lattices of 16 different linear sizes L, from L = 7 to L = 90. 
We studied between 750 samples for L — 90 and 20000 samples for L — 8 and L = 7. 



A. Magnetization 



We study the variation of the absolute value of the magnetization to and of the energy as a function of x. We have 
found that there is a region in x where there are large discontinuities of \m\ and that outside that region to is a smooth 
function of x. The amplitude of the discontinuities is volume independent (this is shown in Fig. ^), while the width 
of the region in x where they appear, shrinks as the volume increases. We analyze our data as follows. For every 
{hi} sample, we chose the rid largest variations of the magnetization between two succesive values of x, Xj = 1/J, and 
Xi + \. Let's call X\ < X2 < ■ ■ ■ < x nd , the values of x at which they occur. The choice of rid is somehow arbitrary. We 
took 2 < rid < 6. It turns out that the x^s fluctuate from sample to sample and their probability distribution is well 
described by a gaussian with mean Xi(L) and variance a(L) 2 . The variance decreases with L and is compatible with 
a power law decay: 

a(L) ~ a L~ s (2) 

We first made the ansatz that: 

Xi = £oo + ~j~ (3) 

i.e. that there is a single discontinuity in the infinite volume limit and that the approach to the thermodynamic limit 
obeys a power law. We obtain a more accurate determination of p by taking differences to eliminate Xao from the fit 

Xij = Xi — Xj — -j— yi) 

We have found that our results cannot be described by a single power for the whole range of sizes we have studied. 
We tried the following alternatives. 

a) Analyze only the largest sizes that are compatible with Eq. ||. We found that sizes larger or equal to 24 3 are well 
described by Eq. ||, with Xoo = 2.26 ± .01, p = .80 and a % 2 per degree of freedom x 2 = 1-29 . With 90% probability 
.76 < p < .83 We also have found the value of the exponent S of Eq. || : .78 < 6 < .86 . The choices of the lower 
size cut-off and of the number of discontinuities rid are quite arbitrary and induce systematic uncertainties on p. In 
order to get an idea on these uncertainties, we tried different "reasonable" choices, i.e. choices with a reasonable 
X 2 per degree of freedom. Keeping only sizes larger or equal to 30 3 , we found that p = .77 and x 2 = -56. With 
90% probability .74 < p < .81. These results have been obtained by analyzing the 5 largest discontinuities of the 
magnetization {rid — 5). We also took rid = 3 and, in this case, .735 < p < .81We conclude that the dependence of p 
on rid is small. 

b) Consider also a subdominant power law correction: 

Xi = Zoo + § (1 + ^ ) (5) 
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This ansatz describes well the data down to a value of L as small as L — 7 with a \ 2 per degree of freedom x 2 = 1. 
and yields p = .66 and q — 1.4 . There are many parameters in this fit, p and q are correlated and the uncertainty on 
p and q much larger than before. The best fit to the data is shown in Fig. ||, while the allowed region in p and q is in 
Fig.|. 

We conclude that the appearance of several discontinuities in the magnetization is a finite volume artifact and that 
our results are fully compatible with the hypothesis of a single discontinuity in the thermodynamic limit. The approch 
to this limit is a very slow power law. The value of this power depends essentially on the assumptions underlying 
the data analysis: a single power law or a power law behaviour with a subdominant correction. This discontinuity 
has not been seen in the previous simulations for the following reason. Only the average magnetization has been 
measured and for finite volumes, the position of the discontinuities fluctuates from sample to sample so that the 
average magnetization seems continuous. 

B. Energy 

Concerning now the energy, it is convenient to separate the energy into two terms, H = —L 3 (JHi + H2), with 

<i,j> i 

For a fixed sample, the value of Hi = E\ (J) in the ground state is the derivative of the energy per spin with respect to 
—J. It can be shown that E\ (J) is a non decreasing function of J, i.e., —H is a convex function of J. We observe large 
discontinuities of Hi as a function of J, as in the case of the magnetization, but, in the present case, the amplitude 
of these discontinuities shrinks as the volume increases. Our data are compatible with the following ansatz: 

A El ~ * (7) 

where AEi, i = 1, •••,?!; denote the ni largest discontinuities of E\. We find a = .48 ± .06 Therefore Hi, i.e. 
the J derivative of the energy, becomes continuous at the thermodynamic limit. Despite the discontinuity of the 
magnetization, the energy and its first derivative are continuous and the exponents are non classical. Eq. ^ sheds some 
light to the long standing question of the existence of different spin configurations almost degenerate in energy. The 
reasonning goes as follows. For J = 0, £1 <~ .0, while for J ^> 1., E\ = 3. In general E\ is a slow varying monotonous 
function of J except in the region .2 < J < .5 where most of the variation occurs. If the spin configuration <r| is a 
ground state for J = with energy per spin — E^ = J^E^ +E^f > and another spin configuration af^ is a ground 
state for J = > J (1 '~> with energy per spin -E^ = J^E {2) + 5 then for J< 3 ' = (E^ - E ( 2 1] ) / (e[ 1] - e{ 2) ), 
the two configurations have exactly the same energy — E^ — J^E± + E^ (J^ < < because of the 

convexity of —H )FL If the largest discontinuity is of the order of L~ a , it follows that in the interval < J < 

(2) (2) 

there are at least ~ (E\ — E\ )L a different configurations with energy E2 < E < E\, which are pairwise exactly 
degenerate for some value of J. This argument is valid for any choice of E\, E2, arbitrarly close to each other. 

C. Disconnected susceptibility 

In order to compare with previous simulations, we have considered the so called disconnected susceptibility Xdis = 
L 3 m 2 i.e. the square of the magnetization averaged over the disorder. It has been assumed that at zero temperature 

Xdis ^ L'-ydx - x^L 1 '") (8) 

where x c is the location of the transition at infinite volume and v the correlation length exponent. In [ p| it was 
found that, for gaussian random fields at zero temperature, fj = 1.1 ± .1 and v = 1 ± .1. In fact Eq. §1 is compatible 



1 One could find the ground state energy for J = j' 3 ', which is either identical to i?*- 3 ' or lower. In the former case the ground 
state energy is exactly known in the entire range J (1) < J < J (2) . In the latter case one could iterate this procedure and so it 
is possible to compute the ground state energy exactly as a function of J. This will be exploited in a forthcoming publication. 
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with our results. To see this, we propose the following simplified model: m = m\ for x < X\ and m = mi for x > x\, 
where x\ is a gaussian random variable with mean x ~ x^ + cjlP and standart deviation a ~ <7qL~ s . A simple 
calculation gives, in this model, Xdis ~ L 3 f((x — Xqq )L S + L s -p). H 8 = p, this agrees with Eq. || provided 77 = 1 
(which is compatible with M]) and <5 = 1/z/. We found above (case a), that 8 = .82± .04 and p = .80± .04, compatible 
with 8 — p. This gives v = 1.22 ± .06 which is slightly different from the value v = 1 ± .1 of [O] . If we restrict our 
analysis to the same (smaller) sizes as in [|llj we find 8 — .96 ± .03, i.e. v = 1.04 ± .03, as in [|11| and p = 1.17 ± .03. 
Alternatively we also have fitted directly our data using Eq. |^. This provides a different method of analysis. Since, 
as already mentionned, a single power is not sufficient to describe the whole range of sizes, we have restricted our 
analysis to the largest sizes. We have found fj = .98 ± .01 in agreement with our simple model and Xoo = 2.265 ± .005, 
a value compatible with the previously deduced value. 

III. BIMODEL DISTRIBUTION OF RANDOM FIELD. 

We consider next the case of a bimodal field distribution, 

P(hi) = (1/2) (8(hi - 1) + 8{hi + 1)) (9) 

i.e. hi — ±1 with equal probability. For every sample {hi}, we chose equally spaced x's with 8x = .0125 and we 
scanned in x. The behaviour of the magnetization is very similar to the case of the gaussian probability distribution 
of the random field. We found again that there is a region in x = hj J where there are large discontinuities of \m\. 
The amplitude of the discontinuities is volume independent, while the width of the region in x where they appear, 
shrinks as the volume increases. We made the same hypothesis as in the case of the gaussian distribution, namely 

a) We found that sizes larger or equal to 24 3 are well described be equ. (3), with x^ = 2.21 ± .01, p = .60 and a x 2 
per degree of freedom \ 2 = -6- With 90% probability .57 < p < .63. Keeping only sizes larger or equal to 30 3 , we 
found that p = .59 and x 2 = -76. With 90% probability .55 < p < .63 These results were obtained by analyzing the 4 
largest discontinuities of the magnetization (rid = 4). 

b) Considering also a subdominant power correction the data are well described with a x 2 per degree of freedom 
X 2 = 1-2 and yields p — .56 and q = 3. Unfortunatly the error bars are much larger than before. The allowed region 
in p and q is presented in Fig. ||. 

For the energy, the situation is more complex. In addition to the type of singularities we observed in the gaussian 
case, there are new stronger singulaties, at integer values of x, i.e. x — 3,4,5, ... These singularities are present for 
all sizes and the amplitude of the discontinuity of E\ at these points seems to be volume independent. Nothing in 
particular is happening to the magnetization at these points. While the presence of these singularities is intuitively 
easy to understand in the case of the bimodal distribution, the decoupling of the magnetization from the energy is 
less intuitive. 



IV. CONCLUSION AND DISCUSSION 

In this paper we have found that he exponent p, which describes the finite size corrections and is usually identified 
with \jv ( v is the correlation length exponent), seems to take different values for the gaussian or the bimodal 
distribution of the random field. This would signal the breaking of universality for the RFIM. It is therefore important 
to discuss a) the uncertainties on the determination of p and b) the identification of p with 1/v. 

a) The statistical errors on p are rather small. The largest uncertainties are due to the assumptions made to analyze the 
data, namely single power law behaviour or inclusion of subdominant corrections in L. But whatever the assumptions 
are, one finds significally different values, provided one uses the same assumptions for both random field distributions. 

b) Let us remind the arguments which lead to the identification of p with 1/V For finite systems of linear sizes L, one 
can define an effective critical temperature T e ff(L). For a second order transition, T e ff(L) could be the temperature 
for which the susceptibility is maximum. The finite size scaling hypothesis asserts that T e //(L) is shifted from T c , 
the critical temperature of the infinite volume system: T e // = T c + cL~ x l v . Assume a second order transition for the 
RFIM and, in agreement with PRG, that the exponents do not depend on the position on the transition line in the 
T J plane, or the direction one crosses it. Assume furthermore, that this is valid down to zero temperature. We then 
expect the locations of the magnetization discontinuities Xi(L)'s to be shifted from their infinite volume value x^: 
Xi(L) = Xac + CiL~ x l v ', i.e. that the exponent p defined above is p — 1/v and therefore is universal. If p is different 
for the gaussian or the bimodal random field distribution as we claimed above, one of the above generally accepted 
assumptions, which can be derived in the context of PRG, in not valid for the random field Ising model. After all it 
is already known that the e expansion leads to incorrect results for the RFIM. The most conservative attitude would 
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be to assume (without any a priori reason) that the bimodal distribution is singular at T = and that the general 
argument does not apply. But we would like to point out that, except from aesthetic appealing, the only theoretical 
argument for universality in disordered systems is based on (unbroken) replicas and perturbative renormalization 
group. And we know that this may be incorrect. We would like to point out that also in the spin-glass case, there 
is numerical evidence that universality may be broken Jlq] ). We feel therefore that universality should be further 
checked for disordered systems. In particular it would be interesting to study other distributions of random fields. 
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FIG. 1. Difference of magnetizations just before xi and just after x$ as a function of 1/L (see text), averaged over random 
field configurations drawn from a gaussian distribution. The error bars are smaller than the symbol. 
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FIG. 3. Figure 3. Allowed regions for the p and q exponents (see text) within a 90% confidence level. The upper left region 
is for a gaussian distribution and the lower right region is for the bimodal distribution of the random fields. 
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